using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;

using Atomic.Libraries;
using Atomic.Libraries.Mathematics;
using Atomic.Libraries.Mathematics.Functions;
using Atomic.Libraries.Mathematics.Optimization;
using Atomic.Libraries.Mathematics.Optimization.Ipopt;

namespace Atomic.Thermodynamics.StateEquations
{
	[Serializable]
	[DebuggerDisplay("{DebuggerDisplay}")]
	public class BirchMurnaghan4 : IStateEquation
	{
		private double e0, v0, b0, bp0;

		public BirchMurnaghan4(double e0, double v0, double b0, double bp0)
		{
			this.e0 = e0;
			this.v0 = v0;
			this.b0 = b0;
			this.bp0 = bp0;

			// The second derivative of the bulk modulus isn't parameter in the BM4 model.
			//bpp0 = (-143.0 + 63.0 * bp0 - 9.0 * bp0 * bp0) / (9.0 * b0);
		}

		public BirchMurnaghan4(Vector e)
		{
			if (e == null)
			{
				throw new ArgumentNullException();
			}

			if (e.Length != 4)
			{
				throw new ArgumentException();
			}

			throw new NotImplementedException();
		}

		public static void Fit(IEnumerable<IVolumeEnergyPoint> points, out double e0, out double v0, out double b0, out double bp0)
		{
			if (points == null)
			{
				throw new ArgumentNullException();
			}

			double bpp0;
			BirchMurnaghan.Solve4(points.ToArray(), out e0, out v0, out b0, out bp0, out bpp0);
		}

		public static FittedBirchMurnaghan4 Fit(IEnumerable<IVolumeEnergyPoint> points)
		{
			double e0, v0, b0, bp0;
			Fit(points, out e0, out v0, out b0, out bp0);

			return new FittedBirchMurnaghan4(e0, v0, b0, bp0, points);
		}

		public static FittedBirchMurnaghan4 Fit(params IVolumeEnergyPoint[] points)
		{
			return Fit((IEnumerable<IVolumeEnergyPoint>)points);
		}
		
		public static FittedBirchMurnaghan4 FitPressure(IEnumerable<IVolumeEnergyPressurePoint> points)
		{
			// Solve standard model first and use parameters as initial values.
			double e0_0, v0_0, b0_0, bp0_0;
			Fit(points, out e0_0, out v0_0, out b0_0, out bp0_0);

			Variable e0 = new Variable();
			Variable v0 = new Variable();
			Variable b0 = new Variable();
			Variable bp0 = new Variable();

			Variable v = new Variable();
			Function e = Energy(e0, v0, b0, bp0, v);
			Function p = -e.Derivative(v);

			Function f = 0.0;
			foreach (IVolumeEnergyPressurePoint point in points)
			{
				f += Function.Sqr(point.Pressure - p.PartialValue(v | point.Volume));
			}

			f = Function.Sqrt(f / points.Count());

			IpoptOptimizer o = new IpoptOptimizer();
			o.Variables.Add(e0, v0, b0, bp0);
			o.ObjectiveFunction = f;
			IOptimizerResult r = o.Run(e0 | e0_0, v0 | v0_0, b0 | v0_0, bp0 | bp0_0);

			return new FittedBirchMurnaghan4(r.OptimalPoint[e0], r.OptimalPoint[v0], r.OptimalPoint[b0], r.OptimalPoint[bp0], points);
		}

		public static FittedBirchMurnaghan4 FitPressure(params IVolumeEnergyPressurePoint[] points)
		{
			return FitPressure((IEnumerable<IVolumeEnergyPressurePoint>)points);
		}

		public static FittedPotentialBirchMurnaghan4 FitPotential(IEnumerable<IVolumePotential> potentials)
		{
			double e0, v0, b0, bp0;
			Fit(FittedPotentialBirchMurnaghan4.CreatePoints(potentials, 0.0), out e0, out v0, out b0, out bp0);

			return new FittedPotentialBirchMurnaghan4(e0, v0, b0, bp0, potentials);
		}

		public static FittedPotentialBirchMurnaghan4 FitPotential(params IVolumePotential[] potentials)
		{
			return FitPotential((IEnumerable<IVolumePotential>)potentials);
		}

		/// <summary>
		/// Energy (eV) as a function of volume (Å^3).
		/// </summary>
		public double Energy(double v)
		{
			// http://en.wikipedia.org/w/index.php?title=Birch%E2%80%93Murnaghan_equation_of_state&oldid=573560674#Expressions_for_the_equation_of_state
			double w = Math.Pow(v0 / v, 2.0 / 3.0);
			return e0 + 9.0 / 16.0 * v0 * b0 * ((w - 1.0) * (w - 1.0) * (w - 1.0) * bp0 + (w - 1.0) * (w - 1.0) * (6.0 - 4.0 * w));
		}

		/// <summary>
		/// Energy as a function of volume.
		/// </summary>
		public Function Energy(Variable v)
		{
			// Same as above, just using abstract (and slower) Function objects.
			Function w = Function.Exp(Function.Log(v0 / v) * 2.0 / 3.0);
			return e0 + 9.0 / 16.0 * v0 * b0 * ((w - 1.0) * (w - 1.0) * (w - 1.0) * bp0 + (w - 1.0) * (w - 1.0) * (6.0 - 4.0 * w));
		}

		public static Function Energy(Function e0, Function v0, Function b0, Function bp0, Function v)
		{
			// Same as above, just using abstract (and slower) Function objects.
			Function w = Function.Exp(Function.Log(v0 / v) * 2.0 / 3.0);
			return e0 + 9.0 / 16.0 * v0 * b0 * ((w - 1.0) * (w - 1.0) * (w - 1.0) * bp0 + (w - 1.0) * (w - 1.0) * (6.0 - 4.0 * w));
		}

		/// <summary>
		/// Pressure as a function of volume.
		/// </summary>
		public double Pressure(double v)
		{
			// http://en.wikipedia.org/w/index.php?title=Birch%E2%80%93Murnaghan_equation_of_state&oldid=573560674#Expressions_for_the_equation_of_state
			//return 3.0 * b0 / 2.0 * (Math.Pow(v0 / v, 7.0 / 3.0) - Math.Pow(v0 / v, 5.0 / 3.0)) * (1.0 + 0.75 * (bp0 - 4.0) * (Math.Pow(v0 / v, 2.0 / 3.0) - 1.0));

			double p, p_v;
			Compute(v, out p, out p_v);

			return p;
		}
		
		public Function Pressure(Variable v)
		{
			return -Energy(v).Derivative(v);
		}	

		public double BulkModulus(double v)
		{
			double p, p_v;
			Compute(v, out p, out p_v);

			return -v * p_v;
		}
		
		public Function BulkModulus(Variable v)
		{
			return -v * Pressure(v).Derivative(v);
		}

		private void Compute(double v, out double p, out double p_v)
		{
			// http://en.wikipedia.org/w/index.php?title=Birch%E2%80%93Murnaghan_equation_of_state&oldid=573560674#Expressions_for_the_equation_of_state

			double w = v0 / v;
			double w_v = -w / v;

			double w23 = Math.Pow(w, 2.0 / 3.0);
			double w23_v = 2.0 / 3.0 * w23 / w * w_v;

			double w53 = w23 * w;
			double w53_v = w23_v * w + w23 * w_v;

			double w73 = w23 * w53;
			double w73_v = w23_v * w53 + w23 * w53_v;

			//p = 3.0 * b0 / 2.0 * (Math.Pow(v0 / v, 7.0 / 3.0) - Math.Pow(v0 / v, 5.0 / 3.0)) * (1.0 + 0.75 * (bp0 - 4.0) * (Math.Pow(v0 / v, 2.0 / 3.0) - 1.0));
			p = 3.0 * b0 / 2.0 * (w73 - w53) * (1.0 + 0.75 * (bp0 - 4.0) * (w23 - 1.0));
			p_v = 3.0 * b0 / 2.0 * ((w73_v - w53_v) * (1.0 + 0.75 * (bp0 - 4.0) * (w23 - 1.0)) + (w73 - w53) * 0.75 * (bp0 - 4.0) * w23_v);
		}

		/// <summary>
		/// Equilibrium energy (eV).
		/// </summary>
		public double EquilibriumEnergy
		{
			get
			{
				return e0;
			}
		}

		/// <summary>
		/// Equilibrium volume (Å^3).
		/// </summary>
		public double EquilibriumVolume
		{
			get
			{
				return v0;
			}
		}

		/// <summary>
		/// Equilibrium bulk modulus (eV/Å^3).
		/// </summary>
		public double EquilibriumBulkModulus
		{
			get
			{
				return b0;
			}
		}

		public double EquilibriumBulkModulusPressureDerivative
		{
			get
			{
				return bp0;
			}
		}

		private string DebuggerDisplay
		{
			get
			{
				return string.Format("E0 = {0}, V0 = {1}, B0 = {2}, B0' = {3}", e0, v0, b0, bp0);
			}
		}
	}
}
